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Abstract 

We perform simulations of an absolute value version of the Villain model on <j> 3 and (f) 4 Feynman 
diagrams, "thin" 3-regular and 4-regular random graphs. The <f> 4 results are in excellent quantitative 
agreement with the exact calculations by Dorey and Kurzepa for an annealed ensemble of thin graphs, 
in spite of simulating only a single graph of each size. We also derive exact results for an annealed 
ensemble of 4> 3 graphs and again find excellent agreement with the numerical data for single 3 graphs. 

The simulations confirm the picture of a mean field vortex transition which is suggested by the 
analytical results. Further simulations on 4> 5 and 6 graphs and of the standard XY model on (f) 3 
graphs confirm the universality of these results. The calculations of Dorey and Kurzepa were based 
on reinterpreting the large orders behaviour of the anharmonic oscillator in a statistical mechanical 
context so we also discuss briefly the interpretation of singularities in the large orders behaviour in 
other models as phase transitions. 



1 Background and Calculations 



It has long been known that mean field behaviour is found in models with short range interactions living 
on tree-like structures such as Bethe lattices Q. This approach circumvents some of the problems that 
appear with using infinite range interactions to get mean field results. Difficulties still arise, however, 
both analytically and numerically when dealing with the dominant boundary of such trees. Random 
graphs, which are locally tree-like and have no external legs, offer a way round this problem, giving a way 
of calculating and simulating models on closed lattices with short range interactions that still behave in 
a mean field like manner. 

There has been a considerable amount of work on spin glass models on random graphs Q mostly 
with Ising or Potts spins, often based on analogy with the corresponding Bethe lattice. Recently it 
was pointed out that transplanting methods from matrix models and 2D quantum gravity allowed a 
considerable simplification of many of the proofs that had been derived || and offered the possibility of 
attacking problems like replica symmetry breaking from a different perspective. In effect one considers a 
d = 0, N = 1 matrix, or scalar, model where N is the size of the matrix, to generate the requisite ensemble 
of random graphs in a Feynman diagram expansion. For large enough graphs saddle point methods can 
be employed and the mean field critical behaviour calculated from the algebraic saddle point equations^]. 

The matrix model inspired approach to discrete spin models was predated by independent work on 
the finite temperature quantum mechanics of the anharmonic oscillator, interpreted as a, d — 1, N = I 
matrix model || . This described a Villain transcription of the continuous spin XY model on 4 graphs 
and gave a mean field vortex transition rather than the Kosterlitz-Thouless (KT) || transition of the 
standard two dimensional XY model. In essence, the thermodynamic limit of the </> 4 random graph model 
is described by large orders in the perturbation series of the anharmonic oscillator in finite temperature 
quantum mechanics. 

Part of the original motivation for considering the XY model transition on an ensemble of random 
graphs was the study of the XY model coupled to 2D quantum gravity. This is equivalent to putting the 
discretized model on an annealed ensemble of fat, or planar, graphs. Calculations JtJ and simulations Q 
have both indicated that the transition was still KT in nature in such an ensemble, but a definitive proof 
of this has been elusive. The idea in Q was to find an ensemble of graphs in which complete calculations 
were possible and in which some sort of transition was still manifest. An ensemble of generic "thin" 
random graphs proved to be such a case. 

If we write the finite temperature partition function for the anharmonic oscillator as 



Z(J3,g) 



D(j>exp — 



dr 



1 



(1) 



where f3 — l/T is the inverse temperature, and carry out the perturbative expansion in g, 

oo 

Z((3,g)^J2 Z k(P)9 k , (2) 



fc=0 



then each Zk(/3) may be written as a sum over Feynman diagrams, which is this case are random <fi 
graphs with k vertices, giving 



z k {p)={-i) k Y J s{G) ... n e ex p 

G ^° ^° i=l <ij>m ij =-oo 



(-\ti - tj + mij(3\) 



(3) 



In the above S(G) is the symmetry factor, which will generically be unity for a large graph, the ti are 
attached to each vertex and integrated over to (3, and the rriij are attached to each link and summed 
over the integers. The partition function Zk(/3) can be thought of as coming from embedding the 4> 
graph on a circle of period (3. The finite temperature one-dimensional propagator that appears in the 
above, 



D <ij>= ^2 exp(-\U-tj +my/3|) 



(4) 



ma=—oo 



1 Inverting the logic of such an approach to use statistical mechanical methods to investigate the large orders behaviour 
of Feynman diagram expansions was actually suggested some years ago by Bender and Wu M . 
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assigns a time coordinate tjj to each end of an edge as well as a winding number my to the edge itself. 
Written in this form the similarity with the Villain version of the XY model || , where the edge factor is 

= £ ex P (-f (*< - 9 J + 27rm y ) 2 ) . ( 5 ) 

niij— — QO ^ ' 

with # taking values in the interval to 2n, is apparent so we might expect the same sort of critical 
behaviour. We have, in effect, interpreted the finite temperature quantum mechanics of the anharmonic 
oscillator as an absolute value version of the Villain model living on thin random graphs. 
For convenience here we perform a rescaling of the i, to obtain 

lk „ „ 

(-/^-tj+myl), (6) 



Z k ((3) = (-l) k (3 k J2S(G) / ... / []^ II E ex P 

G i=l <ij> m i:j =-oo 



which will simplify some of the later formulae for quantities to be measured in the simulations. As in Jj| 
the free energy per vertex is defined as 

F = - Hm llogl^l, (7) 

fe^oo k rik 

where nk is the number of graphs of size k, rik — (16) fe (fc-l)!. Note that, in this case, the thermodynamic 
limit is taken by setting k — > oo by hand. Unlike the case of planar diagrams and two dimensional gravity 
there is no cosmological constant dual to the area that one can tune to induce critical behaviour. 

The energy per vertex in the model p — dF/d/3 measures the expectation value of the target-space 
length of the embedded graph, as can be seen by differentiating the Laplace transform of Zk(0), 

/>oo 

Z fc (/3) = / dLZ k (L)ex.p(-L0), (8) 



with respect to /?. As discussed in || , p can be interpreted as a measure of the density of vortices in the 
model. The specific heat is also given by the standard formula 

d 2 F 

or equivalently by directly differentiating the energy C — dp/dT. 

The analytical solution of the thin graph model proceeds by looking at the large orders behaviour of the 
anharmonic oscillator partition function ]Iq| ] in equ.(l). After a rescaling x = 2^J—g^>, the k — > oo saddle 
point solution is given by the trajectories of period (3 of a particle moving in the potential V = —^x 2 + jX 4 . 
One uses the dispersion relation, 

z M . I /" V mZ( ^ +i0) . do) 

TT J-oo 9 + 

along with the saddle-point evaluation of lmZ(/3,g) to extract Z k {j3). A non-trivial instanton solution 
only exists above a critical value of the period (3 (in the 4> 4 case /3 C = V2ir). For /3 < (3 C the only 
contribution is from the trivial solution x(t) = x m i n , where the particle sits at the bottom of the well. 
The change of behaviour at j3 = j3 c is taken as the signal for a phase transition in the associated absolute 
value Villain model. In the low temperature phase we find 

- (-!) fc ? ( iSH 1/2 inmr ik+1/2) r(* + m, (n) 



2 \2tt dp 



while for (3 < (3 C , 



Z, 
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In the above 1(0) is just the (scaled) classical action of the instanton satisfying 



where E(f3) is the energy associated with this trajectory. To obtain explicit expressions in the above we 
need to determine the dependence of E on the period (3. This can be done by evaluating the first integral 
of the classical equations of motion, 

f X2 dx 

ME) = 2 / - (14) 
J X1 y/2[E-V(x)]' 

where x\ and xi are the turning points. In the </> 4 case, the above integral can be evaluated perturbatively 
near (3 C J5j. Inverting the power series gives E(0) = — 1/4 + 4(/3— /3 c )/(3/3 c ) + . . .. These solutions predict 
p = 1//3 for /3 < (3 C and p ~ exp(— 0) at large (3. Similarly C — 1 for (3 < (3 C and sweeps up to a sharp 
cusp of height 19/3 as (3 — ► /3 C +. In order to obtain the full temperature dependence of p and C shown 
in Figs. 1 and 2, we found it useful to express both f3(E) and the (scaled) classical action in terms of 
elliptic integrals, which can easily be evaluated with high precision. In summary, we see a second order 
transition of the mean-field type, with a sharp cusp discontinuity in the specific heat. 

An explicit calculation can also be carried through in a similar style for the 4> 3 case where the partition 
function is 



Z(fi,g) = J -D^cxp (- J* ' dr 



(15) 



After a rescaling x = —Zgcj), the problem is mapped onto considering the trajectories of a particle in the 



potential = — ix + hx 3 . The relation (13) between the classical action and the energy can be solved, as 



in the </> 4 case, by inverting the relation for (3(E) to get E(j3). Here we find 

/3 = 6 i/ 3 j£2L ) (16) 

where K is the complete elliptic integral with modulus k = (xz — xi)l(xz — x\), and x\ < X2 < x% are 
the roots of 2{E - V(x)\ = 2E + x 2 ~ 2x 3 /3 = 0. Inverting gives E((3) = -1/6 + 6(/3 - /3 c )/(5/? c ) + . . ., 
where (3 C is now 2-7T. The form of p and C is similar to the </> 4 graphs, but we now find a peak of 41/10 in 
C as (3 — * (3 C + and C = 1/2 for f3 < (3 C . As in the 4 case, to compute the full temperature dependence 
of p and C, we first expressed the (scaled) classical action in terms of elliptic integrals, 

1 = i^Vl^ 3 ~ Xl ^' 2 ~ 2k2 + 2k ^ E ^ ~ ( 2 ~ 3fc2 + k4 ) K ( k )] - E P- ( 17 ) 

A couple of other numerical consequences of the analytical results are worth remarking upon as they 
are convenient for verifying that simulations are performing correctly. Firstly, we find that p(T c ) is 
equal to T c on <f> A graphs. Secondly, it should be noted that the (3 C is determined only by the period of 
oscillations around the minimum of the potential for any potential. Approximating this region with a 
quadratic gives (3 C = 2tt/iti if V ~ — Vq + m 2 (x — xq) 2 /2, where Xq is the minimum point. This reproduces 
the explicit results derived above for 3 and <f> graphs with rather less pain. Finally, it is worth noting 
that, for a potential with an anharmonic term of the form 2 (p +1 ) , we will have C = p for T > T c because 
p = pT for T > T c in general. 

We can contrast the critical behaviour described above with the standard XY model on a flat two 
dimensional lattice 

z = n L d ^ ex p ( $ c ° s & - ^) ) ' 

\ <ij> J 





which displays a topologically driven Kosterlitz-Thouless (KT) transition. The specific heat has only a 
broad cusp rather than a divergence. However, the correlation length has a critical singularity 

e = A e ex P (^ F ), (19) 
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as does the spin susceptibility 

X = ^ex P (^ F ). (20) 

The exponent v is predicted to be 1 /2 in the KT theory and the correlation function critical exponent is 
predicted to be 77 = 1/4, where r\ is given by 

x « (21) 

In spite of the differences the physical picture of the transition on thin graphs is still very similar to that 
of the standard KT transition. Given the interpretation of p as a vortex density, the preceding saddle 
point results show that as the temperature is increased (ie (3 is decreased) vortices are liberated, with the 
vortex density increasing by almost an order of magnitude around j3 c . 

As we have already noted, it is not a foregone conclusion that putting the model on any collection of 
random graphs will give mean field behaviour - the model still displays a KT transition on an annealed 
set of planar random graphs (ie when coupled to 2D quantum gravity). An explicit check of the thin 
graph predictions is therefore not totally vacuous. In what follows we describe some modest simulations 
that we carried out in order to verify the mean field nature of the transition for variants of the model on 
various thin graphs. We are essentially interested in the behaviour of the energy and specific heat with f} 
and determining j3 c . We simulated the absolute value Villain action on (j> 3 , (j> 4 , <fi 5 and <fi 6 random graphs, 
as well as looking at the standard XY action of equ.(18) on (j) 3 graphs in order to check it had the same 
critical behaviour. 



2 Simulations 

We need to exercise a little care in defining our observables in the simulation because of the unusual form 
of the Boltzmann factors in the partition function pT| . If we define the auxiliary sums 

00 

£ = 2J exp - t o + m v\) ' 

7rtij— — OQ 

00 

Ei = ^2 \U - tj + m i: j I exp (-/?|t< - tj + my]) , (22) 

rriij — — 00 
00 

S 2 = \U - tj + rriijf exp (-f3\U - t 3 + m lj \) , 

iriij — — 00 

then the definition of the energy p — dF/df3 applied to the partition function in equ.(6) gives 

for the energy per site, where the < > denote a thermal average and the additional 1 / (3 comes from the 
overall factor of (3 k that appears when the U are rescaled. 

The specific heat can be obtained in the simulations either by direct numerical differentiation of the 
measurements of p using C — dp/dT, or by differentiating equ.(6) twice to give 

C = ^(</>-< p >V^(e(!-!))-1. (24) 

The second term is non-canonical and is due to the summed Boltzmann factors whereas the —1 results 
from the overall (3 k . 

Having decided on our observables, it now remains to choose an update scheme for the simulation. 
A simple Metropolis update can be used quite efficiently by adopting a discrete approximation to the 
periodic tj and then tabling the resulting Boltzmann factors and associated sums So, 1,2, which can then 
be looked up during the course of the simulation [fi"l|| . Depending on the temperature, we took from 100 
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to 1000 different t values for the tables, which were constructed by truncating the sum over rriij at ±100. 
Increasing these limits made no appreciable difference to the measured results in any of the simulations 
reported here. In doing this we are taking a Zioo • • • ^1000 approximation to the 0(2) symmetry of the 
model. Notice that for the absolute value version of the Villain model one has to be quite careful with this 
approximation, in particular at low temperatures, since the associated discretization error enters linearly 
in the action and not squared as in the standard Villain model. It is perhaps worth remarking that one 
could equally well envisage leaving the m,j as free variables on each link to be sampled in the course of 
the simulation, but previous work in which the sum is carried out a priori as here has given good results 



for the standard Villain/XY model |11 and we stick to this. The direct form of D<y> converges rapidly 



for large /3, but is also possible to use a "dual" representation that is best suited for small (3 simulations, 

^•> = |fl+E^^^|. (25) 
P \ k y =i l-i < 2Trfe 



The numerical differences, even when using the representations well into the "wrong" regions of (3 for 
both, are slight. 

The final ingredient in the simulations is the choice of a random graph. The calculations we have 
outlined in the first section are supposedly for an annealed ensemble of thin graphs, so in theory we should 
carry out "flip" moves in the same fashion as in simulations of 2D gravity on planar graphs. However, we 
can evade this responsibility by appealing to previous simulations of the Ising model on random graphs 
H , where a single random graph of a given size was enough to extract the mean field critical behaviour 
predicted in an annealed calculation. As the thin random graphs are essentially tree-like it appears that 
one tree is very much like another as far as a ferromagnetic transition is concerned, although a (quenched) 
average over graphs becomes essential when non-self-averaging transitions such as those in spin glasses 
are simulated. We shall take the lazy man's approach here and verify that the single graph results are 
consistent with the analytical calculations that, strictly speaking, only apply to an annealed ensemble. 

We simulated the absolute value Villain model of equ.(6) on <j) 3 and </> 4 graphs of size up to 2500 
vertices, as well as small runs of N — 250 (f> 5 and <fi 6 graphs. In all the cases we carried out 500,000 
Metropolis sweeps at each (3 value, with a measurement every tenth sweep, after allowing a suitable 
amount of equilibration time. As we have indicated no flip moves were carried out on the graphs con- 
cerned. The energy and specific heat, defined as above, were the principal observables. In addition, we 
also simulated the standard XY model of equ.(18) on 3 graphs, using a single cluster update, largely as 
a check on the universality of the Villain model results. 

Turning now to the results themselves we can see in Fig. 1 that the energy matches closely the 
predictions of |^] . There is a low temperature exponential growth with a "knee" at the phase transition 
T = 1/V27T (~ 0.225) followed by linear growth at larger T. It is also clear from Fig. 1 that p(T c ) = T c , 
as predicted by the analysis in the previous section. We can obtain the specific heat either from direct 
differentiation of the energy or by measuring the observable defined in equ.(24), and both are plotted in 
Fig. 2 for a graph of size 1000. This is the smallest graph size at which the peak reaches its expected 
value (19/3) - larger graphs give similar results, whereas the peak is appreciable lower and more rounded 
on the smaller graphs simulated. Surprisingly, it is the numerical differentiation that gives the closest fit 
to the analytical results, perhaps an indication that our statistics around the phase transition point are 
rather too modest for complete numerical accuracy in quantities involving variances, such as the directly 
measured specific heat. 

Fortified by the good agreement between the analytical and numerical results for 4 graphs we can 
move on to look at the possible variations on the theme that were outlined in the introduction. We 
consider the (f> 3 graphs first. In Fig. 3 we plot the energy for various graph sizes, showing clearly the 
similarity with the (j) 4 results. The "knee" in the curve is at the expected value of T c — 1/2tt ~ 0.16 g. As 
our data points are rather scarcer than for the 4> A graphs, the specific heat peak obtained by numerical 
differentiation as shown in Fig. 4 is not so convincing, but it is clearly of the same general form and has 
the correct large T limit of 1/2. 

We have not carried out such extensive simulations of the cf> 5 and <fi 6 random graphs, simply contenting 
ourselves with verifying that the general form of the energy is similar and that the large T limit is correct. 

2 In this case, however, we do not expect that p(T c ) = T c , but rather p(T c ) = T c /2. 
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In Fig. 5 the energy is plotted up to very large T for 4> 3 ,<fi 4 , <fi 5 and 4> 6 graphs of size 250. From the slopes 
it is clear that the specific heat prediction C — > p for T — > oo on 2 (p +1 ) graphs is satisfied to a high 
degree of accuracy. The finite size effects for a given graph size increase with the degree of the vertices, 
which agrees with the intuitive picture of 2 ( p+1 ' graphs being more "tree-like" for smaller p with a given 
number of vertices. 

All of the simulations and discussions so far have dealt with the absolute value Villain model on various 
random graphs. We have done this because the edge factor in this case is just the one dimensional finite 
temperature propagator of equ.(4), which allows us to borrow (steal!) the results from large orders 
expansion of various anharmonic oscillators. We would expect that the standard XY model of equ.(14) 
would still give us similar results on grounds of universality, but as we have no analytical calculations 
to fall back on in this case it is worthwhile verifying this explicitly with simulations. We therefore 
simulated the standard XY model on cj) 3 graphs of various sizes, with similar statistics to the Villain 
model simulations but using a single cluster update for improved efficiency. 

The specific heat, measured directly in the simulation, is plotted in Fig. 6 for various lattice sizes, 
where it is clear that, although the small and large T limits are different from the Villain models (1/2 
and respectively), there is still a sharp cusp in the curve. This would again indicate a transition of 
mean field rather than KT type, where there is a much gentler bump in the specific heat curve away 
from the phase transition point. 

As for the standard Villain model |jl^| , the differences can be understood by an approximate mapping 
of the absolute value version of the Villain model onto the cosine model. By adapting the formulas in 
Ufl to the present case, we find that the temperature scales should be related by Ji(/3 cos )//o(/3 cos ) = 
(P/2tt) 2 / (l + (/3/27r) 2 ), where 7o,i are Bessel functions and /3 cos denotes the inverse temperature of the 
cosine model. Inserting (3 C = 2n this predicts h(^ os ) / I ((3 C C OS ) = 1/2 or T c cos = 0.8625..., in good 
agreement with the peak location observed in Fig. 6. 



3 Conclusions and Other Models 



The saddle point predictions for the energy and specific heat of the absolute value Villain model on various 
random graphs are verified by the simulations we have carried out. As one might expect the standard 
XY model on "thin" tfi 3 graphs behaves in an analogous fashion, with a mean-field-like transition rather 
than a KT transition. The second point worth emphasizing is that we have not needed to simulate an 
annealed ensemble of random graphs to get good quantitative agreement with the theory, just as for the 
Ising ferromagnet on thin graphs This should be contrasted with the planar graphs in 2D gravity 
where an annealed sum, usually implemented by flip moves in a simulation, appears to be essential. It 
appears that, unless one is very unlucky, any large thin graph is as good as another. Trees, and the near 
trees that thin random graphs represent, all really do look much alike. 

The methods used in deriving the analytical results are borrowings from standard large orders lore 
in quantum mechanics. It is rather the interpretation in terms of a statistical mechanical model that is 
novel, just as the behaviour of discrete spin models on random graphs is extracted from known saddle 
points results for ordinary integrals in [||. We are not restricted to the simple anharmonic oscillator 
in searching for statistical mechanical interpretations of large orders behaviour in quantum mechanics. 
Another possible example was considered by one of the authors of this paper in [ p"3| , namely the quantum 
mechanics of an anisotropic anharmonic oscillator, where the partition function is 



Z 



Dcf>iD(j)2 exp 



dr 



1 



\\ + 4>l + 4>\ + <j>\) + 9(4 + 2c4>\4>l + <t>i) 



(26) 



Taking c as the control parameter rather than 0, the large orders behaviour shows a transition at c = 1, 
where the model is rotationally symmetric. For — 1 < c < 1 the quartic 
instanton solution is 4>\(t) = u(t),4>2{t) = 0, where 



b\ terms dominate and the 



u(t) 



1 



2|A| cosh(i-i ) : 



(27) 



whereas for c > 1 the <\>\<\)\ term dominates and the solution is of the form (f>i{t) = (t) — u(t)/\/2. 
Looking at the Feynman diagrams generated by the model we can see that it is a sort of loop gas, with (f>i 
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loops and <p2 loops mixing via the <j>\<j>2 vertex and a propagator D < ij > between the individual vertices 
on loops of both types. A matrix model style order parameter can be defined as 



M = frf (28) 

which, using the instanton solutions, gives M = 1 for — 1 < c < 1 and M = for c > 1. Alternatively 
we could consider the ratio of mixed to pure vertices, 

M — I— 2 ^ (29) 

with similar results. With either order parameter we can see that the random graphs are filled with x 
loops only for —1 < c < 1 and an even mixture of x and y loops for c > 1. The singular behaviour at 
c = 1 can thus be viewed as a sort of magnetization transition. 

It would be an interesting exercise to see if other large orders results in quantum mechanics could be 
recast in a statistical mechanical mould as in this case. We leave this for future work. 
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Figure Captions 

Fig. 1. The energy for various cj) 4 graph sizes. A dotted line indicates the analytical prediction. 

Fig. 2. The specific heat for a <j> 4 graph of size 1000, obtained via both numerical differentiation of p 
and direct measurement. Again the dotted line represents the analytical curve. 

Fig. 3. The energy for various <fi 3 graph sizes. The dotted line indicates the analytical prediction. 

Fig. 4. The specific heat for cj) 3 graphs of various sizes, obtained via numerical differentiation of p. The 
dotted line shows the analytical prediction. 

Fig. 5. The energy for (f> 3 , cf> 4 , (f> 5 and 4> 6 graphs. The linear prediction, p = pT on </> 2 ( p+1 ' graphs when 
T >T C , is shown as dotted lines to emphasize the very good fit. 

Fig. 6. The specific heat for the standard XY model on <f) 3 graphs of various sizes. 
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